Interacting Anisotropic Dirac Fermions in Strained Graphene and Related Systems 
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We study the role of long-range electron-electron interactions in a system of two-dimensional 
anisotropic Dirac fermions, which naturally appear in uniaxially strained graphene, graphene in 
external potentials, some strongly anisotropic topological insulators, and engineered anisotropic 
graphene structures. We find that while for small interactions and anisotropy the system restores the 
conventional isotropic Dirac liquid behavior, strong enough anisotropy can lead to the formation of a 
quasi-one dimensional electronic phase with dominant charge order (anisotropic excitonic insulator) . 
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Since the isolation of graphene [lj, a two dimensional 
(2D) allotrope of carbon, there has been substantial the- 
oretical and experimental effort to understand and uti- 
lize its remarkable mechanical , thermal Q , electronic 
[3] and transport [f| properties. The charge carriers in 
a pristine graphene sheet move on a honeycomb lattice 
(Fig. Ufa)) as if they were massless relativistic parti- 
cles. The physical properties at low energies are gov- 
erned by the Dirac Hamiltonian resulting in a linear dis- 
persion and isotropic cones (with circular cross section) 
near the Dirac points (Fig. H£b)). Mechanical defor- 
mations, such as uniaxial strain for example along the 
armchair direction (Fig. [He)), can lead to anisotropic 
dispersion with the formation of elliptical Dirac cones 
(Fig. [TJd)) [6]. Within the tight-binding and ab ini- 
tio schemes, the anisotropic Dirac dispersion is found 
for weak to moderate uniaxial strain in any direction, 
while for large strain the electronic structure becomes 
very different along the armchair and zig-zag directions, 
ultimately leading to cone merger and gap formation in 
the latter case [f| . It has also been shown that anisotropic 
Dirac cones can be formed by applying an external peri- 
odic potential [?J on graphene. Such changes of the band 
structure can provide exciting possibilities for "strain en- 
gineering," i.e. manipulation of graphene's electronic, 
optical, etc. properties by applying lattice deformations 
or potentials Q. Strongly anisotropic Dirac cones can 
also appear in certain topological insulators Q, and, in 
a recent development, highly tunable honeycomb optical 
lattices and molecular graphene systems have been cre- 
ated providing a possible route towards exploring 
various anisotropic phases. 

The subject of the present work is the interplay be- 
tween Dirac fermion anisotropy and electron interactions. 
It is known that near the Dirac point, close to charge 
neutrality, unscreened long-range electron-electron inter- 
actions in isotropic graphene can manifest themselves in 
a variety of ways [111 ]: perhaps most spectacularly inter- 
actions lead to logarithmic renormalization of the elec- 
tron spectrum leading to reshaping of the Dirac cones 
(velocity increase), observed in undeformed suspended 
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graphene 12]. Previously anisotropic Dirac fermions 
have been studied in QED3-type models, relevant for the 
cuprate superconductors, where the fermion anisotropy 
was found to be irrelevant in renormalization group (RG) 
13(, i.e. the systems flows towards the isotropic 
This behavior can be reversed, i.e. the anisotropy 
increases, if the fermions couple to a nematic order pa- 
rameter, as suggested for the cuprates fli| . 

In this work we show that in graphene-based 
anisotropic Dirac systems, as they arise in modified (i.e. 
strained, or artificially engineered) graphene, a rich va- 
riety of behavior can take place due to electron-electron 
interactions. Within the RG approach we find that for 
small interactions and anisotropy, the system flows to- 
wards the conventional isotropic fixed point. However 
the anisotropy also favors a transition towards an ex- 
citonic insulator (whose exsistance is well established 
in the isotropic case [Hj]), making this phase accessi- 
ble for much smaller interactions compared to isotropic 



(a) 





- a i 



























(b) 





FIG. 1: (Color online) Honeycomb lattice for (a) undeformed 
and (c) unixially deformed graphene along the y (armchair) 
direction. The corresponding isotropic (b) and anisotropic (d) 
Dirac cones are shown. 
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graphene. We show that provided the anisotropy is large 
enough, a quasi one-dimensional phase with dominant 
charge density wave order (anisotropic excitonic insula- 
tor) can form even for small interactions. Thus strongly 
anisotropic graphene-based systems can relatively easily 
experience interaction-driven transitions towards novel 
electronic phases, providing a promising avenue for ex- 
ploration of unconventional many-body physics. 

Dispersion anisotropy. We consider, for dehnitiveness, 
graphene under uniaxial strain. Without interactions, 
the Hamiltonian is 0: Ho — v x p x a x + v y p y a y , where 
a x , <7 y are the usual (2 x 2) Pauli matrices. Here v x = 
t2a x y / Arf^l, v y = t2a y are the velocities along the x— 
and y— direction respectively; r) = |^ = |j, a x — a/2 and 

a y = -^a, with a ps \/3xlA2 A, as shown in Fig. [Ha). 
The electronic dispersion, as shown in Fig. QJd) is: 



with N being the number of fermion flavors (N — 4 for 
graphene.) Using the standard decomposition E*- 2 ) = 
£(7o£o + ''/ /'/ .;■ + VyPy&yTiy we obtain the dressed 
Green's function 



G(p,e) = 



Z(l) 



SO-Q - V X {l)p x a X - Vy{l)PyGy 



(3) 



where v x (l)/v x = Z(l)(l + E B (0), v y {l)/v y = Z(l)(l + 
Ej,(i)), and the quasiparticle residue is Z(l) = (1 — 
Eo(O)" 1 wl + So(0- At low energies, |p| = p -» 0, all 
quantities diverge logarithmically in terms of the param- 
eter I — ln(A/p), where A is the ultraviolet cutoff. We 
obtain by direct calculation the divergent contributions: 



E o (0 = -al 



N 
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In(A/p), 



(4) 



e(p) = ± 



(i) 



and we have also defined, and will use from now on the 
anisotropy parameter 8. We will explore the range — 1 < 
8 < 0, i.e. < v y /v x < 1, which in the notation of 
Fig.QJc) would imply strain in the armchair (y) direction, 
with v y < v x . Clearly 8 = corresponds to the isotropic 
case, while 8 = — 1 is the limit of decoupled chains. Strain 
in the zig-zag (x) direction can be easily described by the 
same parameter range with appropriate relabeling of the 
axes. The anisotropy parameter 8 is proportional to the 
strain From now on we use the effective anisotropic 
dispersion {T]) which, as mentioned earlier, is relevant to a 
variety of graphene related systems, and is not necessarily 
due to strain. 

Electron- electron interactions and anisotropy. The 
bare long-range Coulomb potential in graphene is given 
by V(p) = 7rF§p where k is the appropriate dielectric 
constant. We consider graphene at charge neutrality 
where the chemical potential \x = and V(jp) is un- 
screened. The interaction effects can be incorporated 
into the self-energy E(p,e), so that the fermion Green's 
function (GF) is given by: G _1 (p,e) = e&o — v x p x a x — 
VyPy&y — E(p,e), where &q is the identity matrix. We 
work in the two-loop approximation, i.e. E(p,e) = 

£ (2) (p,£) = *E q /_°li7^ 2) (p-q,£-^G( )(q, w ) 

with the effective interaction V^(q, uj) = V(q, u>) + 
(U(q,w)) 2 n(q,w), and G (0) is the free GF. The physi- 
cal reasons behind the use of the two-loop approxima- 
tion will be discussed later and vertex corrections are 
neglected as their effect is small, similarly to the case 
of isotropic graphene fl6| . The dynamical polarization 
bubble for anisotropic Dirac fermions is easily evaluated 
to be: 



E*,v(0 = {^*Ii,3 - 2A{ i + 6 f l ( 3 - 4I 2 , 4 )|ln(A/p). 

(5) 

Here Ij = h(8), i = 1,2,3,4 are evaluated as follows. De- 
fine C(9, 8) = cos 2 6» + (1 + 8) 2 sm 2 6, then we have I n (S) = 
i r 27T cos 2 0de f i o „„j t rx\ — i r 27r sin 2 ede 



2tt Jo [C(0,<5)] n/2 ' 2tt JO ^/ci^S)' 

T m _ 1 r 2 ^(l+S) 2 sin 2 6de 

^\°) — 2^ Jo c(eJ) • 

Since 1 + 8 = v y jv x < 1 , it is convenient to define the 
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dimensionless Coulomb interaction coupling a x = 
From Eqs. (I4I5[) we obtain the RG equations for a x (l) 
and the anisotropy 8(1): 



for n = 1, 2, and l 3 (8) = ± J Q 



| = (1 +8)^(h(8) -h(S))- 



-^(1 2 (S)-I 4 (8)) (7) 



We assume 1 + 8 > 0. It is also instructive to write these 
equations in the small anisotropy limit, 8 <C 1: 



rf7 



In the opposite limit of strong anisotropy, <5 ~ — 1, we 
obtain: 



da,, a 2 2a| ^ a 

■dT = ~V + — ' H = {1 + s) - ln \i + 8 



3 

(10) 



n(q,w) 



A r 



2 2 



16t> x v„ 



w 2 g 2 



« 2 (? 2 



RG results: possible electronic phases. The numerical 
(2) solution of the RG equations, Eq. and (JTJ), with N = 
4, leads us to the RG flow in Fig. [2] which is the main 
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FIG. 2: (Color online) RG flow in the plane (a x ,S), in terms 
of the interaction a x — and the anisotropy, 5 = Vy v — . 
There is a stable fixed point at (0, 0) and unstable fixed points 
at (0, — I) and (3/4,0). The broken line a x ,c(5) separates a 
region of a x flowing to zero, for a x < a XtC (8), and the region 
ot x > &x,c(S) where a x either increases or flows towards a 
finite value. The evolution of the couplings at the four points 
labeled by numbers are further shown in Fig. [3] 



result of this work. The arrows represent the variations 
of the couplings as the RG parameter I — ln(A/p) varies 
from towards oo (low-energy limit). 

First, notice that at zero anisotropy (5 — 0) there is 
an unstable fixed point at a c — 3/4; the flow towards 
strong coupling at a x > a c corresponds to an excitonic 



insulator [15[ , i.e. a system with a finite gap in the spec- 
trum (see later discussion.) Within the perturbative RG 
scheme this fixed point arises from the competition be- 
tween the one and two loop contributions [16[ , as is clear 
from Eqs. (JSJ) , d9j) - Of course one cannot hope to reliably 
obtain the exact value of a c at intermediate coupling, 
but the very existence of a c has been well established by 
different methods 15J, with the result a c « 1.1. Thus 
we use the two-loop result in a x and then take into ac- 
count the anisotropy exactly. The dashed line a XlC (8) in 
Fig. [5] separates a phase with a flow towards the stable 
fixed point at a x — 0, 5 = and a region with a flow to- 
wards large (diverging) or finite a x . We will show shortly 
that the latter case, a x > a x ^ c (S), is indeed characterized 
by a divergent susceptibility towards excitonic insulator. 
For small S we find that a x (l) becomes large beyond the 
modified excitonic transition line, which can be calcu- 
lated from Eq. © to be: a x<c (6) « (3/4)(l-|5|/4). Thus 
the anisotropy shifts the excitonic transition to smaller 
coupling (a XfC (5) decreases as \S\ increases). 

Even below the excitonic transition line a x < ct XjC (5) 



the anisotropy flow of 6(1) towards zero can be very slow 
and non-monotonic, as shown in Fig. [3] This is especially 
true for initial values close to the excitonic line. Notice 
that 8(1) in fact increases for small / up to I = 10 — 
20. Therefore such a signature can be easily observed for 
graphene at finite Fermi energy (or temperature). For 
example if at T — the Fermi energy is [i ~ 0.2 meV 
we have to stop the RG flow at I* = \n(v x A/n) ps 10 
(taking v x A — 5 eV), resulting in an increase of 6 up to a 
factor of 2 (or more) depending on the proximity to the 
excitonic line (curves 1,2 in Fig. [3]). 
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FIG. 3: (Color online) Evolution of the couplings a x (l) (upper 
panel) and S(l) (lower panel) for initial values corresponding 
to the four numbered points from Fig. [2] 

For a x > a x . c (S) the region of moderate anisotropy 
and relatively weak coupling, 0.4 < \S\, a x < 0.6, be- 
low the dashed line (Fig. [2]) deserves special attention. 
Here the flow is towards finite a x and 5 — > — 1, i.e. the 
one-dimensional (ID) limit. From Eq. (fTD|) . valid close 
to the ID line, we can see that a* = 3/(27r) = 0.48 
separates regimes where a x either experiences (relatively 
small) decrease, or increases under RG flow as the sys- 
tem approaches the ID limit. Flow of \S\ beyond 1 is not 
allowed since it reverses the sign of the velocity. 

In order to characterize the nature of the phase for 
&x > ct x>c (6), we introduce in the Hamiltonian infinites- 
imal couplings A M to potential order parameters (in the 
particle- hole channel) : A M J^k ^k^M^k, M = X, y, z. We 
then compute the interaction corrections to the corre- 
sponding vertex functions T M = A^(l + m (A/p)) and 
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find the susceptibilities: 



J Z (S), Xx.y = Xz 



L -W*)- (11) 



r 2 *' cos 2 ede 



Jx(5) — 2ir Jo [C(6>,5)] 3 / 2 

and C(9, 5) = cos 2 + (1 - 



Here J,(S) = * J^*L=, 

J v\°! — 2tt Jo [C(e,5)] 3 / 2 ' 
(5) 2 sin 2 #. The quantities Xm(0 depend on the scale I 
via the RG running of a x (l),S(l). Divergence of Xm(0 
at I = I* signals spontaneous breakdown of symme- 
try in the corresponding channel (l7j . We have found 
that Xz(0 i g th e fastest divergent susceptibility through- 
out the region a x > a x ^ c (5), and therefore the sys- 
tem is an excitonic insulator with an order parameter 
($r^a z i&) = (i> A 'ipA) — {WBi } B) 7^ which describes charge 
density modulation between the two sub-lattices . In 
Fig. H] we plot Xz(0 within the perturbatively accessible 
parameter regime when the system flows towards the ID 
limit, and a x is finite. In this case Xz diverges logarith- 
mically: 

Xz(l) ~ MO/*] ln(4/[l + 6(1)}), 1 + 5(1) -> 0. (12) 

One can estimate the transition temperature T c from the 
formula I* — ln(w x A/T c ), which gives for example (see 
Fig.HJl: T c w 2K (I* w 10.3), T c w 45K (Z* w 7), T c « 5x 
10 3 K (/* 2.2). Naturally T c increases with increasing 
a x . 



x 




FIG. 4: (Color online) Dominant divergent susceptibility 
Xz(l) corresponding to the order parameter a^iS) , for dif- 
ferent initial couplings (a x > a X:C (8)). The divergence is at 
1 = 1* where the system approaches the ID limit, 5(1) — 1. 



weak coupling while non-perturbative methods (such as 
excitonic pairing equations as well as bosonization) are 
needed to further quantify the exact shape of the phase 
boundary in the strong coupling regime (a x ~ 1). How- 
ever the overall topology of the phase boundary as well as 
the quantitative behavior at small a x are well captured 
within our approach. 

Finally, in the weak coupling region a x < 0.5 we find 
that the renormalization factor Z decreases as the system 
approaches ID. From Eq. (|4j we obtain ^ = — gn+|) i 
which leads to the following behavior near the ID line: 



Z(l) ~ (1 + <5(Z)) 1/4 , 1+*(0->0 



(13) 



Notice that Z vanishes with a universal exponent 1/4, 
independent of the interaction a x . 

Summary and outlook. In conclusion, we have found 
that the interplay between the anisotropy of the Dirac 
spectrum and long-range electron-electron interactions 
can lead to rich variety of behavior and two main elec- 
tronic regimes: (i) a weak coupling phase character- 
ized by a flow toward isotropic Dirac physics but with 
strongly renormalized parameters at low energies, and 
(ii) a quasi one-dimensional phase with dominant charge 
density wave order (anisotropic excitonic insulator) . The 
strong anisotropic tendencies can manifest themselves 
in a variety of ways, such as squeezing of the Landau 
level spectrum in magnetic field [l9j], increase of the 
density of states which would favor itinerant ferromag- 
netism and affect (increase) the specific heat [ll|, and 
strong anisotropies in transport 6]. The value of the 
interaction on typical substrates is a ~ 0.5 — 0.9 [llj ]. 
while strain of about 15% in the zig-zag direction leads 
to anisotropy \5\ ps 1/2 Q, sufficient to push graphene 
into the anisotropic regime. We also expect that artifi- 
cially engineered graphene (optical lattices and molecular 
graphene [n| ) can provide a promising way to investigate 
the strongly anisotropic states found in the present work. 
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From Eq. (|llj) we find numerically that Xy also diverges 
in the ID limit, but is always smaller than Xz, mean- 
ing that nematic-type (gapless) order is competing in the 
ground state. It is interesting to note that in studies of 
sliding Luttinger liquid (SLL) phases [l8| one typically 
finds smectic non-Fermi liquid metals as well as charge 
ordered states. Notice that the irrelevance of the inter- 
chain hopping (our v y ) is one of the main characteristics 
of the SLL. Our analysis is valid, strictly speaking, for 
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